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, We perform quantum Monte Carlo simulations to study the 2D hard-core Bose-Hubbard 



model in a random potential. Our motivation is to investigate the effects of randomness on 
the Kosterlitz-Thouless (KT) transition. The chemical potential is assumed to be random, 
by site, with a Gaussian distribution. The KT transition is confirmed by a finite-size analysis 
of the superfluid density and the power-law decay of the correlation function. By varying the 
variance of the Gaussian distribution, we find that the transition temperature decreases as 
C . the variance increases. We obtain the phase diagram showing the superfluid and disordered 

"j^ , phases, and estimate the quantum critical point (QCP). Our results on the ground state 

C/3 , reveal the existence of the Bose glass phase. Finally, we discuss what the value of the variance 

, at the QCP indicates from the viewpoint of percolation. 

KEYWORDS: Kosterlitz-Thouless transition, superfluidity, quantum Monte Carlo simulation, 
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o 

■ 1. Introduction 

> ' 

ff^ ' The Kosterlitz-Thouless (KT) transition^) is an important concept underlying two- 

. dimensional (2D) physics. The XY model is a typical model describing the KT transition. 

cn 



This model describes well not only 2D XY magnets and 2D crystals but also superfluidity in 
thin films, which is a particularly good example to use when studying the KT transition. One 



o 

' of the characteristics of the KT transition can be seen in superfluid (SF) density. Kosterlitz 



- ^ . and Thouless theoretically predicted that the SF density ps shows a discontinuous jump at 

XI ■ 

' the transition temperature Tkt given by the normalized form 



2 

Ps = -Tkt- (1) 

vr 

Equation (1), which is known as the universal jump, provides an important proof of the KT 
transition. Many experiments using a torsional oscillator^' have observed this universal jump. 
Another characteristic is the correlation function of the order parameters. According to the 
KT theory, the phase transition results from unbinding vortex pairs. When the temperature T 
is lower than the transition temperature Tkt ■, the correlation function of the order parameters 
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shows a power-law decay as a function of the distance r, {'ip{x)ip{x + r))oc r^''>^'^\ which is 
called the quasi long-range order (QLRO). This state is stable because the excitation by 
vortex pairs does not destroy the correlation of the order parameters. As T increases above 
Tkt, vortex pairs become unbound as free vortices. These free vortices destroy the QLRO, 
forcing the system into the disordered phase, where the correlation decays exponentially. 

Finite-size effects must be taken into account in numerical simulations of the KT tran- 
sition. Although finite-size scaling is the best analysis method for studying thermodynamic 
limits with finite-size numerical data, ordinary finite-size scaling cannot be applied to numer- 
ical simulations of the KT transition^^ because the logarithmic correction of the system size 
makes it difficult to obtain reliable estimates of Tkt- Several numerical studies^' have used 
the XY model to analyze the nature of the KT transition. These previous studies focused 
on the SF density to obtain Tkt • Weber and Minnhagen performed classical Monte Carlo 
simulations for the classical 2D XY model and used the KT renormalization group equation to 
analyze the numerical data instead of using ordinary finite-size analysis.^) They systematically 
studied the size dependence of the SF density at the transition temperature to reveal from the 
logarithmic dependence that the KT transition actually occurs in the 2D classical XY model. 
Harada and Kawashima verified the KT transition of the 2D quantum XY model using quan- 
tum Monte Carlo (QMC) simulations. They also calculated the SF density precisely and 
confirmed a close agreement between their numerical data and the logarithmic scaling of the 
KT theory. 

For about two decades, the properties of disordered Bose systems have been an interesting 
topic of theoretical and experimental study. The Bose-Hubbard model gives a good description 
of the physics of interacting bosons, 

<i,j> i i 

where S|(Si) is a boson creation (annihilation) operator, = S|Si, J is the hopping param- 
eter, U is the on-site repulsive interaction, and fXi is the random chemical potential at site 
i. This model corresponds to the XY model when we impose the hard-core Hmit {U ^ oo). 

The hopping parameter, the interaction of particles, and the disorder of the chemical poten- 
tial compete to determine the state of Bose systems. After the theoretical prediction of Bose 
glass (BG),'^-' which is the intermediate state between the SF state and the insulator state, 
theoretical^"^ and experimental research^^^^^^ has been done to study disordered Bose sys- 
tems extensively. Experimentally, such systems can be realized by confining liquid ^He into 
nano-porous materials, such as Vycor,-*^^) Aerogel,^^) and Gelsil.^^^ These systems can be also 
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realized by trapping atomic gases in an optical lattice. The optical lattice is a powerful 
stage on which to study lattice Bose systems because the parameters described in eq. (2) can 
be easily controlled. 

This paper describes numerical simulations of 2D disordered Bose systems. Although 
there are many studies of these systems focused on the ground state/°' -"^^^ they include little 
discussion on the relationship between the KT transition and randomness. Our motivation is 
to investigate this relationship while taking proper account of finite-size effects and how the 
KT transition can survive against randomness. The following is an outline of this paper. In 
Section 2 , we describe the model and the numerical method. We perform QMC simulations to 
study the disordered BH model for a 2D square lattice with a periodic boundary condition. The 
chemical potential differs from site to site according to a Gaussian distribution. In Section 
3, calculations and results are shown. The KT transition is confirmed from the finite-size 
scaling of the SF density and the power-law decay of the correlation function. To investigate 
the effects of randomness, we vary the variance of the distribution. Although the transition 
temperature decreases as the variance of the distribution increases, the KT transition still 
survives in the presence of randomness. The quantum critical point (QCP) at which the SF 
phase disappears at T = is estimated. By calculating the SF density and the compressibility, 
the nature of the ground state in this disordered system is investigated. Our results reveal 
the existence of the BG phase. In Section 4, we discuss several implications of the estimated 
QCP value. Our results can be understood from the viewpoint of percolation. Finally, Section 
5 is the Conclusion . 

2. Model and numerical method 

We study the KT transition in lattice Bose systems using eq. (2) with the hard-core limit 
(U-^oo). This limit implies that the number hi of particles at site i is restricted to or 1. 
Then, eq. (2) reduces to 

<i,j> i 

This model is equivalent to the 5=1/2 quantum XY model in a random field. We apply eq. (3) 
to a 2D square lattice of size Ns=I? with a periodic boundary condition. In this study, the 
random chemical potential fii is assumed to obey a Gaussian distribution given by 

where /xq is the mean and is the variance . In our calculations, J is fixed at 0.5 by the 
correspondence with the S=l/2 quantum XY model,^) and the mean /xq is —1.0. The effects 
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of randomness are investigated by varying the variance o"^. 

The computational method used is path-integral QMC simulations using a modified 
directed-loop algorithm, which is one of the most efficient for calculating the lattice BH 
model. We perform calculations for a maximum system size L = 250 and maximum 
number of random sampling 512. 

3. Calculation and results 

In this section, we calculate the SF density ps, the two-body correlation functions, and 
the compressibility k. The main purpose is to investigate whether the KT transition occurs 
even in these disordered systems. The analysis procedure is the following. We perform the 
finite-size scaling of ps to determine precisely Tkt- Then, we assess whether the correlation 
function below Tkt shows a power-law decay. The variance is varied from 0.0 to 7.0 to 
study the effects of randomness on the KT transition. We use finite-size analysis in the same 
way as in the case of the pure KT transition in ref.5,6. Even if the system is disordered, the 
symmetry of eq. (3) should be maintained, so that we can assume the finite-size analysis for 
the pure KT transition remains valid. In Section 3.4, the temperature dependence of the SF 
density and compressibility are shown to reveal the existence of the BG phase. 

3. 1 SF density and estimation of the transition temperature 

The SF density ps is essential for studying the KT transition because it shows the universal 
jump at Tkt- It can be calculated directly by estimating the winding number W using path- 
integral QMC simulations. The SF density ps is defined as 

Ps = ^{W'). (5) 

As an example. Fig. 1(a) shows the temperature dependence of the SF density Ps{L) for 
various system sizes L = 24, 32, 48, 64, 96, and 128 with = 1.5, allowing us also to confirm 
the size dependence of the SF density. When the system size is relatively small, the SF density 
increases gradually as temperature decreases. As the system size increases, the SF density rises 
abruptly at a certain temperature. This seems to support the universal jump at Tkt in this 
disordered system, if we extrapolate to the case in which L is taken to infinity. 

What we want to know is the transition temperature Tkt- The solid line in Fig. 1(a) 
represents ps = (2/7r)T/J. The intersection between the solid line and /Os(oo) = lim.L^^ps{L) 
gives Tkt for the thermodynamic limit. The problem is estimating Ps{oo) from data on Ps{L). 
As mentioned above, the results of numerical simulations are affected by finite-size effects. 
In order to obtain a reliable Tkt, we follow the method of former studies^' to perform 
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the finite-size analysis of Ps{L) for this disordered system. According to these studies, the 
finite-size correction of at T = Tkt can be given by 

p,(L) = g^(l+ / ), (6) 

vr V 21og(L/Lo(rKT))^ 

where Lq is the characteristic length of the order of the lattice constant. It is notoriously 
difficult to determine Tkt precisely from this finite-size correction;^) however, we propose a 
simple method to obtain Tkt precisely by plotting 1/{psTt/T — 2) as a function of L. 
We rewrite eq. (6) as 

logL - logLo(rKT) = -JTT-L ^- (7) 

According to eq. (7), we assume that l/(ps(L)-7r/rKT — 2) should behave like log L— log Lo(Tkt) 
as a function of L. Figure 1(b) is a rescaled plot of the SF density of Fig. 1(a) showing 
l/(/9s7r/r — 2) as a function of L for various temperatures. The bold curve in Fig. 1(b) 
represents \ogL — log Lq{T]<:y) ■ Here, we regard Lq as a fitting parameter because changes of 
Lq only shift the plot of log L— log i>o(7kT) versus L vertically without changing its shape. We 
vary Lq to obtain the value at which the plot of log L — log Lo(^kt) fits the rescaled data of 
1/{PsTt/T-2). When Lq is 0.90(±0.02), the plot of log L-logZo(rKT) fits T/J = 0.166(±0.02) 
(See Fig. 1(b)). Thus, the transition temperature Tkt with = 1.5 is estimated to be 
0.166(ib0.02) J. This logarithmic dependence of Ps{L) at exactly Tkt is the key to confirming 
the KT transition. We perform this finite-size scaling to determine Tkt(<7^) for cr^ from 0.0 
to 7.0 and confirm the logarithmic behavior of the KT transition. 

3. 2 Correlation functions 

An important characteristic of the KT transition is the existence of the QLRO showing 
the power-law decay below Tkt- We calculate the two-body correlation function of a distance 
rto investigate the power- law below Tkt, as follows: 

G{r) = ^Y.^h^Mr). (8) 
^ i 

Figure 2 shows the correlation functions with cr^ = 6.0 for three temperatures, namely 
T=Tkt(~0.0268 J), T > Tkt (0.5 J), and T < Tkt (0.0067 J), for a finite system size (T = 128). 
The critical exponent of the correlation functions r]{TKT) is determined to be ^ theoretically 
and numerically. It is clear that the numerical data supports r]{TKT)=\ from the solid line 
in Fig. 2 and the power-law decay below Tkt- For each value of cr^, the correlation function 
below Tkt decays according to the power law. This indicates the existence of the QLRO and 
also confirms the existence of the KT transition even in disordered systems. 



5/13 



J. Phys. Soc. Jpn. 



Full Paper 




Fig. 1. (a) (color online) Temperature dependence of the superfluid density for various system sizes 
with CT^ = 1.5. The size dependence of the SF density is apparent. The lines connecting points are 
guides to the eye. The solid line is {2/tt)T/ J. (b) (color online) Rescaled plot of Fig. 1(a). The 
sohd curve is the fit to logL — logLo(7kT)- 




Fig. 2. (color online) Correlation functions with — 6.0 for various temperatures. We can observe 
the power-law decay of the correlation functions below Tkt~0.0268 J. 



3.3 Phase diagram and estimation of the QCP 

Figure 3(a) shows the phase diagram in the a^-Tj J plane. The points correspond to the 
numerical data of TxtI^^) obtained by QMC simulations, and the solid line connecting the 



6/13 



J. Phys. Soc. Jpn. Full Paper 

points is a guide to the eye. The transition temperature Tkt decreases as increases, but 
the SF phase survives within a broad region of cr^ . This indicates that the SF phase is robust 
against randomness. 

This robustness can be understood from the viewpoint of the Harris criterion. The 
Harris criterion is apphed when considering the relevancy of randomness to the second-order 
phase transition. Near the critical temperature T^, the specific heat Cobeys C oc |Tc/(T— Tc)|", 
and the critical exponent a is a key to the Harris criterion. According to the Harris criterion, 
randomness is relevant when the critical exponent a is positive. The properties of disordered 
systems qualitatively differ from those of pure systems. On the other hand, the randomness 
is irrelevant when the critical exponent a is negative (the specific heat does not diverge); the 
criticality does not change even in disordered systems. In the present case, since we study 
the KT transition, the specific heat does not diverge and shows a broad peak. Although the 
critical exponent a for the KT transition cannot be defined, randomness can be expected to 
be irrelevant to the KT transition if we consider the non-diverging specific heat. 




I ' ■ ■ ■ ' ' ' ■ 1 
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Fig. 3. (color online) (a) Phase diagram in the ct^- Tj J plane. The solid line is a guide to the eye. (b) 
A rescaled plot of Fig. 3 (a). The quantum critical point (cr^) is estimated as «8.79. 

Here, we define the QCP where the SF phase disappears. To estimate the critical variance 
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(Jc^, we rcplot the data of Fig. 3(a) in the T/J-(l/cj^) plane (Fig. 3(b)). The data in the 
low-temperature region can be fitted to (l/o"^ — l/alY, and the fitting parameters, l/o"^ and 
X, are estimated to be 0.114(±0.003) and 0.946(=b0.04), respectively. Thus, the value of is 
estimated to be 8.79(±0.03). 



3.4 Properties of the system for > u\ 



The nature of the quantum phase transition in the disordered Bose system has been a 
major topic of debate. As shown in rcf.7, in the groimd state, there are three phases in the 
disordered Bose system: the SF phase, the insulator phase, and the BG phase. According 
to ref.7, in the pure system, the phase transition between the SF and the insulator phases 
occurs directly. However, in the disordered system, there is no direct transition between the 
SF and the insulator phases because of the intermediate state, namely the BG state. The 
physical variables necessary to distinguish these three phases are the SF density and the 
compressibility. The SF (insulating) phase is characterized by finite (zero) SF density and 
compressibility. Although the insulating phase has vanishing compressibility, the insulating 
BG phase has finite compressibility. 

In the previous section, we found the QCP and that the value of a\ is 8.79. In order to study 
the properties of the system for > we calculate the SF density and the compressibility. 
Figure 4(a) shows the temperature dependence of the SF density with = 0.0, 2.0, 4.0, 
and 15.0 (L = 40). As increases, the SF density decreases. Clearly, the SF density with 
(T^ = 15.0 vanishes. 

The compressibility k is defined as 

1 

Here, N is the total number of particles. Figure 4(b) shows the temperature dependence of 
the compressibility. The data shown in Fig. 4(b) are for = 0.0, 2.0, and 15.0, /xq = —1.0, 
and L = 40. In addition, we plot the data of the insulating phase with cr^ = 0.0, /xq = 
—5.0, and L = 40. The compressibility is calculated for various system sizes, but little size 
dependency of compressibility is apparent. In the insulating phase, for the lower temperature, 
the compressibility vanishes. As the temperature increases, the compressibility has an thermal- 
activated form e~'^ao.vlT' ^ith a finite energy gap Agap-^^-* With \xq = —1.0, when is 0.0, we 
can see that the compressibility has a sharp peak. The peak becomes broad as increases, 
but the compressibility remains finite. With = 15.0, although the system seems to be in 
the insulating phase because of the zero SF density, the compressibility remains finite even if 
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the temperature becomes zero, just as in the case <t^ < o"^. This reveals the existence of the 
BG phase in this disordered system. 





Fig. 4. (color online) (a) Temperature dependence of the SF density for various variances with L = 40. 
(b) Temperature dependence of compressibility with L — 40. 



4. Discussion 

In this section, we consider what the computed value of cr^ indicates. We consider the 
phase diagram of the ground state in a pure system ((T^=0.0 and jjLi=^ in eqs. (2) and (3)) 
analytically and numerically to study the relation between the SF phase and the chemical 
potential. This relation leads us to compare the present problem with the 2D site-percolation 
problem. 

Analytical approach: In our calculations, the hard-core limit was imposed on eq. (2). This 
limit implies that the number of particles at each site is limited to or 1. The phase boundary 
between the insulator and the SF phases for the pure system can be obtained easily under this 
restriction. For large ^, the ground state is the insulator phase given by \'4jq) = |1, 1, . . . , 1) = 
Y\f=i in which a single particle occupies each site. The first excited state given by iV'i) = 
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Z^fcLi Wk) is the superposition of |afc), where |ofc) is the state in which a site k is vacant 
while other sites {i^ k) are occupied. Comparing the energy between the ground and the first 
excited state, we can obtain the phase boundary of = 4 J in the 2D system. 

Numerical approach: We investigate the chemical potential dependence of the average of 
the number of particles at each site n=^. Figure 5 shows an example for hopping parameter 
J = 0.5 at various temperatures. As mentioned above, the transition point /Xc between the 
insulator and the SF phases is defined as the point at which n changes from unity. The 
temperature dependence of fic is apparent from the figure. By extrapolating ndT) as a function 
T, //c(0) can be determined. 



0.995 
0.99 
0.985 
0.98 
0.975 

-2.5 -2.4 -2.3 -2.2 -2.1 -2 

\^ 

Fig. 5. (color online) Chemical potential dependence of the average number of particles at each site 
for J — 0.5 at various temperatures. 
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In this way, we obtain the phase diagram of the ground state in the pure system (Fig. 6). 
The dotted lines are the analytical results = 4J, and the points are the numerical results. 
Clearly there is a nice agreement between these results. The central region is the SF phase, 
and the upper and lower regions are the insulator phases. 

In our study, the SF state can be described as the state that the correlation of particles 
develops over the whole system, namely, a particle moves over the whole system during the 
imaginary time (3 development. We consider what the value of a1 indicates from the viewpoint 
of percolation by comparing Fig. 6 and QMC simulations for our disordered case at T = 0. In 
QMC simulations for the disordered case, the random distribution of the chemical potential 
was assumed to be Gaussian distribution. The mean was —1.0 and the variance cr^ was 
calculated to be 8.79. The hopping parameter J was fixed at 0.5. As shown in Fig. 6, for 
the pure system, the system shows superfluidity in the region of -2< /i <2 at J = 0.5. The 
chemical potential /x = ±2 can be considered as the critical chemical potential for systems 
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Fig. 6. (color online) Phase diagram of the ground state in a 2D pure system. The dotted lines are 
the analytical results and the points are the numerical results. The central region is the superfluid 
phase, and the upper and lower regions are the insulator phases. 



to show superfluidity. Next, the probability Pq that the chemical potential at a site takes 
a value within -2< /i <2 when the random distribution of the chemical potential obeys the 
Gaussian distribution {fiQ = —1.0, o"^ = 8.79) is calculated to be 47.6%. If the probability 
that the chemical potential takes a value within -2< fi <2 exceeds Pq, the system can enter 
the SF phase. This probability Pq should be considered as the threshold for entering the SF 
phase. 

The values of probability Pq reminds us of the 2D site-percolation problem. It is the 
geometrical connection that is important for the site-percolation problem. In this picture, the 
sites exist randomly. That is, we assign only two energy levels, zero or infinity, randomly to 
the site energy. When the probability that the zero-energy sites are distributed exceeds the 
2D site-percolation threshold Pc, which was obtained to be 59.3% numerically, a spanning 
cluster exists and the system becomes coherent. If we map our disordered case to this two- 
level system, sites such that -2< fi <2 might be considered as the zero-energy sites and the 
others as the infinite-energy sites. However, no precise correspondence between Pq and Pc 
can be seen. Apparently, Pq is about 10% smaller than Pc- Therefore, we wish to know what 
is underlying this difference. 

There are two reasons why Pq and Pc should be different. The first reason is the difference 
of the percolation process. If the percolation process is regarded as dynamical, meaning that 
particles propagate infinitely far from their initial position, whether the dynamics of the par- 
ticle is classical or quantum makes a clear difference to the percolation process. The threshold 
Pc = 59.3% was obtained under the classical picture. In this picture, the geometrical connec- 
tion can be understood as the propagation of the particles through zero-potential sites. On the 
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other hand, in the quantum pereolation deseription, not only the site energy but also the state 
of the site, namely, the number of particles hi for each site, must be considered. Our numerical 
results for the disordered systems refer to the limit of the dynamical percolation process along 
imaginary time as /3(=l/T)^oo, and the state of the site is the number of particles. Thus, 
our case should be dominated by the quantum percolation process. The second reason is that 
the site energy is not infinite but has a finite value in the present case. Under the quantum 
picture, the infinite site energy blocks propagation of the quantum particle, causing a quan- 
tum localization effect. Because of this localization effect, the quantum percolation threshold 
is larger than Pc in this case.^^^ However, it is known that the finite site potential reduces the 
quantum percolation threshold below Pc- If the site energy is finite, the quantum particle can 
tunnel through the potential barriers, and it tends to propagate infinitely far from the initial 
position, unlike in the classical description.^^) In our case, the site energy takes a finite value 
that obeys a Gaussian distribution. Together, these two reasons explain why Pq should be 
expected to be smaller than Pc- 

5. Conclusion 

In this paper, we described, using QMC simulations, a study of the 2D hard-core BH model 
with a random chemical potential that was assumed to obey a Gaussian distribution. The KT 
transition was confirmed from a finite-size analysis of the SF density ps and the power-law 
decay of the correlation function. In order to investigate the effects of randomness on the 
KT transition, the variance of the distribution was varied, and we obtained a phase diagram 
showing the SF and disordered phases. Even if the variance increases, the KT transition 
survives robustly against the randomness. The Harris criterion helps to explain this robustness 
of the KT transition. We observed the QCP and estimated the value of to be 8.79. The 
nature of the phase for a'^>af. revealed the existence of the BG phase. 

In addition, we considered what this value of cr^ indicates and argue that it can be un- 
derstood from the viewpoint of percolation. Comparing the 2D percolation threshold and the 
present case makes the quantum effects in the percolation process clear. Our results should 
be dominated by the quantum percolation processes. 
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